function [PSDPeaks, fPeaks] = psdpeaks(data, fS, nMajorPeaks, freqRes, minPeakDist, minPeakHeightRatio, doPadding, doDetrend, doPlot)
%PSDPEAKS 计算信号PSD的主要峰值
%
% Input Arguments
% # data: 向量，信号
% # fS: 标量，采样频率，单位Hz
% # nMajorPeaks: 标量，返回的主要峰值个数，默认为5
% # freqRes: 标量，频率分辨率，单位Hz，默认为0.01Hz
% # minPeakDist: 标量，峰值之间的最小间隔，默认为1Hz
% # minPeakHeightRatio: 标量，峰值距最小值的高度（相对于峰峰值的比例），默认为0.2
% # doPadding: 逻辑标量，true表示当所有峰值数量不足nMajorPeaks时在返回值后添加NaN以使其长度等于nMajorPeaks，默认为true
% # doDetrend: 逻辑标量，true表示去除data中的常数项，默认为false
% # doPlot: 逻辑标量，true表示绘图，默认为true
%
% Output Arguments
% # PSDPeaks: 长度为nMajorPeaks或峰值总个数（当doPadding为false且峰值总个数小于nMajorPeaks时）的列向量，峰值的PSD值，单位dB/Hz
% # fPeaks: 长度与PSDPeaks相同的列向量，峰值的频率值，单位Hz

if nargin < 3
    nMajorPeaks = 5;
end
if nargin < 4
    freqRes = 0.01;
end
if nargin < 5
    minPeakDist = 1;
end
if nargin < 6
    minPeakHeightRatio = 0.2;
end
if nargin < 7
    doPadding = true;
end
if nargin < 8
    doDetrend = false;
end
if nargin < 9
    doPlot = true;
end

if doDetrend
    data = detrend(data, 'constant');
end

nWnd = 2^nextpow2(fS/2/freqRes);
[PSD, fPSD] = pwelch(data, nWnd, [], nWnd, fS);
PSDLog = log10(PSD);
[~, peakLocs] = findpeaks(PSDLog, 'MINPEAKDISTANCE', round(minPeakDist/(fS/2/length(fPSD))), ...
    'MINPEAKHEIGHT', min(PSDLog)+minPeakHeightRatio*peak2peak(PSDLog), 'SORTSTR', 'descend');
nActualMajPeaks = min(length(peakLocs), nMajorPeaks);
if (nActualMajPeaks<nMajorPeaks) && doPadding
    PSDPeaks = NaN(nMajorPeaks, 1);
    fPeaks = NaN(nMajorPeaks, 1);
end
PSDPeaks(1:nActualMajPeaks, 1) = PSD(peakLocs(1:nActualMajPeaks));
fPeaks(1:nActualMajPeaks, 1) = fPSD(peakLocs(1:nActualMajPeaks));

if doPlot
    semilogy(fPSD, PSD); % 绘制PSD图
    hold on;
    plot(fPSD(peakLocs), PSD(peakLocs), 'k^', 'markerfacecolor', [1 0 0]); % 绘制PSD峰值
    for iPeaks = 1:nActualMajPeaks % 标注主要峰值
        text(fPeaks(iPeaks), PSDPeaks(iPeaks), [num2str(fPeaks(iPeaks), '%.2f'), ', ', num2str(PSDPeaks(iPeaks), '%.2f')]);
    end
    grid on;
    xlabel('Frequency(Hz)');
    ylabel('PSD(dB/Hz)');
    hold off;
end
end